#include <iostream>
#include <vector>
#include <cmath>
#include <matplot/matplot.h> // 需安装matplot++库用于绘图

using namespace std;
using namespace matplot;

// 生成连续正弦信号：y = sin(2πft)，f为信号频率
vector<double> generate_continuous_signal(double f, double duration, double dt) {
    vector<double> t, y;
    for (double t_val = 0; t_val < duration; t_val += dt) {
        t.push_back(t_val);
        y.push_back(sin(2 * M_PI * f * t_val));
    }
    return y;
}

// 对连续信号采样
vector<pair<double, double>> sample_signal(const vector<double>& t_cont, const vector<double>& y_cont, double f_s) {
    vector<pair<double, double>> samples;
    double dt_sample = 1.0 / f_s; // 采样间隔
    double t_prev = -dt_sample;
    for (size_t i = 0; i < t_cont.size(); ++i) {
        if (t_cont[i] - t_prev >= dt_sample - 1e-6) { // 按采样间隔取点
            samples.emplace_back(t_cont[i], y_cont[i]);
            t_prev = t_cont[i];
        }
    }
    return samples;
}

// 用线性插值重建信号（简化重建，实际可用低通滤波）
vector<double> reconstruct_signal(const vector<pair<double, double>>& samples, const vector<double>& t_cont) {
    vector<double> y_recon;
    size_t idx = 0; // 采样点索引
    for (double t : t_cont) {
        // 找到当前t所在的两个采样点之间
        while (idx < samples.size() - 1 && samples[idx + 1].first < t) {
            idx++;
        }
        if (idx >= samples.size() - 1) {
            y_recon.push_back(samples.back().second);
            continue;
        }
        // 线性插值
        double t0 = samples[idx].first, y0 = samples[idx].second;
        double t1 = samples[idx + 1].first, y1 = samples[idx + 1].second;
        double y = y0 + (y1 - y0) * (t - t0) / (t1 - t0);
        y_recon.push_back(y);
    }
    return y_recon;
}

int main() {
    // 信号参数：频率f=10Hz（最高频率），时长1秒，连续信号时间步长0.001秒（高频采样，近似连续）
    double f_signal = 10.0;
    double duration = 1.0;
    double dt_cont = 0.001;
    vector<double> t_cont;
    for (double t = 0; t < duration; t += dt_cont) {
        t_cont.push_back(t);
    }
    vector<double> y_cont = generate_continuous_signal(f_signal, duration, dt_cont);

    // 情况1：采样频率f_s1=25Hz（≥2*f_signal=20Hz）
    double f_s1 = 25.0;
    auto samples1 = sample_signal(t_cont, y_cont, f_s1);
    vector<double> y_recon1 = reconstruct_signal(samples1, t_cont);

    // 情况2：采样频率f_s2=15Hz（<2*f_signal=20Hz）
    double f_s2 = 15.0;
    auto samples2 = sample_signal(t_cont, y_cont, f_s2);
    vector<double> y_recon2 = reconstruct_signal(samples2, t_cont);

    // 绘图对比
    subplot(2, 1, 1);
    plot(t_cont, y_cont, "b-", "Original Signal");
    hold(on);
    vector<double> t_samples1, y_samples1;
    for (auto& p : samples1) { t_samples1.push_back(p.first); y_samples1.push_back(p.second); }
    plot(t_samples1, y_samples1, "ro", "Samples (f_s=25Hz)");
    plot(t_cont, y_recon1, "g--", "Reconstructed");
    title("Sampling at f_s ≥ 2f (Nyquist Condition)");
    xlabel("Time (s)"); ylabel("Amplitude");
    legend();

    subplot(2, 1, 2);
    plot(t_cont, y_cont, "b-", "Original Signal");
    hold(on);
    vector<double> t_samples2, y_samples2;
    for (auto& p : samples2) { t_samples2.push_back(p.first); y_samples2.push_back(p.second); }
    plot(t_samples2, y_samples2, "ro", "Samples (f_s=15Hz)");
    plot(t_cont, y_recon2, "g--", "Reconstructed (Aliased)");
    title("Sampling at f_s < 2f (Aliasing Occurs)");
    xlabel("Time (s)"); ylabel("Amplitude");
    legend();

    show();
    return 0;
}